library(velocyto.R)
library(SeuratWrappers)
library(Seurat)
source("tianfengRwrappers.R")
preprocess <- function(loom_path, n){
velodata <- ReadVelocity(file = loom_path)
#匹配两次的barcode
tfunc <- function(s)
{
s <- strsplit(s,".*:",fixed = F)[[1]][2]
s <- strsplit(s,"x",fixed = T)[[1]]
s <- paste(s,as.character(n),sep = "-")
return(s)
}
velodata[["spliced"]]@Dimnames[[2]] = as.character(lapply(velodata[["spliced"]]@Dimnames[[2]],tfunc))
velodata[["unspliced"]]@Dimnames[[2]] = as.character(lapply(velodata[["unspliced"]]@Dimnames[[2]],tfunc))
velodata[["ambiguous"]]@Dimnames[[2]] = as.character(lapply(velodata[["ambiguous"]]@Dimnames[[2]],tfunc))
return(velodata)
}
ctrl <- readRDS("ctrl.rds")
ctrl_velo <- as.Seurat(x = preprocess("./HSD0524.loom",1))
ctrl_velo <- subset(ctrl_velo, cells = WhichCells(ctrl))
ctrl_velo <- ctrl_velo %>%
PercentageFeatureSet(pattern = "^MT-", col.name = "percent.mt") %>%
SCTransform(vars.to.regress = "percent.mt", verbose = F,assay = "spliced") %>%
RunPCA() %>% FindNeighbors(dims = 1:20) %>%
RunUMAP(dims = 1:20) %>%
FindClusters(resolution = 0.1)
# for(i in levels(Idents(ctrl))){
# (DimPlot(ctrl_velo, reduction = "umap",cells.highlight = WhichCells(ctrl, expression = `celltype` == i), pt.size = 0.5, label = T)+ggtitle(i)) %>% print()
# }
ctrl_velo@reductions[["umap"]] <- ctrl@reductions[["umap"]]
ctrl_velo@reductions[["umap"]]@cell.embeddings <-
ctrl@reductions[["umap"]]@cell.embeddings[colnames(ctrl_velo),]
DimPlot(ctrl_velo, reduction = "umap",group.by = 'seurat_clusters', pt.size = 0.5, label = T)
ctrl_velo$Classification1 <- ctrl$celltype
Idents(ctrl_velo) <- ctrl_velo$Classification1
ctrl_velo@assays[["RNA"]] <- ctrl@assays[["SCT"]]
# saveRDS(ctrl_velo,"ctrl_velo.rds")
# ctrl_velo$Classification1 <- Idents(ctrl)
ident.colors <- colors_list[1:length(x = levels(x = ctrl_velo))]
names(x = ident.colors) <- levels(x = ctrl_velo)
cell.colors <- ident.colors[Idents(object = ctrl_velo)]
names(x = cell.colors) <- colnames(x = ctrl_velo)
png("velocity.png")
show.velocity.on.embedding.cor(emb = Embeddings(object = ctrl_velo, reduction = "umap"), vel = Tool(object = ctrl_velo, slot = "RunVelocity"), n = 200, scale = "sqrt", cell.colors = ac(x = cell.colors, alpha = 0.5),
cex = 0.8, arrow.scale = 3, show.grid.flow = TRUE, min.grid.cell.mass = 0.5, grid.n = 40, arrow.lwd = 1,
do.par = FALSE, cell.border.alpha = 0.1)
dev.off()
h5ad_output <- "hsd_velo.h5ad"
library(reticulate)
library(Matrix)
writeMM(t(ctrl_velo@assays$RNA@counts), file='combined.mtx')
writeMM(t(ctrl_velo@assays$spliced@counts), file='spliced.mtx')
writeMM(t(ctrl_velo@assays$unspliced@counts), file='unspliced.mtx')
write.csv(rownames(ctrl_velo@assays$spliced@counts), file = "genes.csv", row.names = FALSE)
write.csv(ctrl_velo@reductions$umap@cell.embeddings, file = "umap.csv", row.names = FALSE)
write.csv(ctrl_velo@reductions$pca@cell.embeddings, file = "pca.csv", row.names = FALSE)
write.csv(colnames(ctrl_velo@assays$spliced@counts), file = "cells.csv", row.names = FALSE)
write.csv(ctrl_velo@meta.data, file = "meta.csv", row.names = FALSE)
source_python('~/scRNAseq/vascular-analysis/build.py')
build(h5ad_output, pca = TRUE, umap = TRUE)
filename = "hsd_velo.h5ad"
import anndata
import pandas as pd
from itertools import chain
#import scvelo as scv
ad = anndata.read_mtx('spliced.mtx')
ad.layers['spliced'] = anndata.read_mtx('spliced.mtx').X
ad.layers['unspliced'] = anndata.read_mtx('unspliced.mtx').X
ad.obs = pd.read_csv('meta.csv',header=0)
ad.obsm['X_umap'] = pd.read_csv('umap.csv',header=0).values
ad.obsm['X_pca'] = pd.read_csv('pca.csv',header=0).values
ad.var_names = pd.Index(list(chain.from_iterable(pd.read_csv('genes.csv', header=0,dtype='string').values.tolist()))).astype(str)
ad.obs_names = list(chain.from_iterable(pd.read_csv('cells.csv', header=0).values.tolist()))
ad.write(filename)
file.remove('combined.mtx')
[1] TRUE
file.remove('spliced.mtx')
[1] TRUE
file.remove('unspliced.mtx')
[1] TRUE
file.remove('genes.csv')
[1] TRUE
file.remove('cells.csv')
[1] TRUE
file.remove('umap.csv')
[1] TRUE
file.remove('pca.csv')
[1] TRUE
file.remove('meta.csv')
[1] TRUE
Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).
The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.